function fp = jacob(f,x,f0)
% fp = jacob(f,x,fo)
%    this function produces a crude forward finite difference
%    approximation for the function f(x) about x; fo is the initial
%    evaluation of f
n   = length(x);
%f0  = feval(f,x);
eps = 1e-8;
fp  = zeros(n);
for j=1:n
    v           = zeros(n,1); 
    v(j,1)      = 1.0;
    xx          = x+v*eps;
    f1          = feval(f,xx);
    fp(:,j)     = (f1 - f0)/eps;
end